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Abstract 

A numerical procedure has been developed to investigate the nonlinear and strain rate dependent 
deformation response of polymer matrix composite laminated plates under high strain rate impact 
loadings. A recently developed strength of materials based micromechanics model, incorporating a set of 
nonlinear, strain rate dependent constitutive equations for the polymer matrix, is extended to account for 
the transverse shear effects during impact. Four different assumptions of transverse shear deformation are 
investigated in order to improve the developed strain rate dependent micromechanics model. The 
validities of these assumptions are investigated using numerical and theoretical approaches. A method to 
determine through the thickness strain and transverse Poisson's ratio of the composite is developed. The 
revised micromechanics model is then implemented into a higher order laminated plate theory which is 
modified to include the effects of inelastic strains. Parametric studies are conducted to investigate the 
mechanical response of composite plates under high strain rate loadings. Results show the transverse 
shear stresses cannot be neglected in the impact problem. A significant level of strain rate dependency 
and material nonlinearity is found in the deformation response of representative composite specimens. 


Introduction 

There is a growing need in military and civil applications for composite materials that not only have 
good structural characteristics, but also good penetration resistance and greater strength after impact. 
Polymer matrix composites are very susceptible to projectile impact such as fragments, flying debris, or a 
failed rotor blade. If not contained, a projectile traveling at ballistic velocity could penetrate a polymer 
matrix composite plate or shell and cause fires, hull damage, occupant injury and component malfunction. 
For example, an engine containment system must be capable of containing fragments traveling at ballistic 
velocities from a failed fan blade. Also, such an impact event may compromise the structural integrity of 
the composite and lead to catastrophic failure. Thus, it is important to develop the ability to predict the 
deformation and failure behavior of polymer matrix composites subject to high strain rate loading 
conditions. An analysis methodology must be able to account for any strain rate effects, material 
nonlinearities, and transverse shear effects that may be present in the impact problem. Also, the 
computational efficiency is critical in the numerical analysis of such a problem. The objective of this 
work is to develop an efficient numerical framework for the analysis of engine containment systems made 
of polymer matrix composites addressing the issues discussed above. 
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Several equivalent single -layer laminate theories have been developed to investigate the mechanical 
response of composite laminated plates on a macroscopic scale. 1 As is well known, Classical Laminated 
Plate Theory (CLT) neglects the transverse shear stresses by using the Kirchhoff hypothesis. By relaxing 
the Kirchhoff hypothesis, i.e., by allowing the transverse normals to not necessarily remain perpendicular 
to the deformed midsurface of the composite laminate, the First-Order Shear Deformation Laminated 
Plate Theory (FSDT) was developed. 2 In FSDT, transverse shear strains are represented as constant 
through the laminate thickness and shear correction factors are used to model the transverse shear 
deformation. By assuming a cubic displacement field, a more accurate theory, Higher Order Laminated 
Plate Theory (HOT), was developed. 3 4 In HOT, transverse shear stresses are modeled by assuming cubic 
displacement fields through the laminate thickness, thereby eliminating the need for shear correction 
factors. Also, unlike FSDT, the conditions of zero transverse shear stresses on the top and bottom surfaces 
are satisfied by imposing these criteria as constraints. 

Various micromechanics models have been developed to characterize the nonlinear properties of 
polymer matrix composites. In micromechanics, the overall properties and response of the composite are 
computed based on the properties and response of the individual constituents. Plasticity and 
viscoplasticity based constitutive equations can be used to compute the nonlinear response of the matrix 
constituent (assuming the fiber is linear elastic, as is usually the case for polymer matrix composites), and 
homogenization methods can then be applied to compute the overall (nonlinear) response of the 
composite. The methodology is based on defining a unit cell in the composite, for which the behavior of 
the unit cell can be assumed to be equivalent to the response of the entire composite ply. Several 
examples of this approach have been developed previously. For example, the Method of Cells (MOC) 
was developed by Aboudi for unidirectional fiber-reinforced composites whose constituents were elasto- 
viscoplastic materials. 5 Pindera and Bednarcyk reformulated the MOC using simplified uniform stress 
and strain assumptions, which resulted in improving the computational efficiency of the analysis 
procedure considerably.' 1 By applying a similar approach and discretizing the composite unit cell into 
three subcells, a two dimensional model that described the elastic -plastic behavior of fibrous composites 
was developed by Sun and Chen. 7 This model then was extended to three dimensions by Robertson and 
Mall. 8 A more precise elastic micromechanics model was proposed by Whitney in which the unit cell was 
divided into an arbitrary number of rectangular, horizontal slices. 9 Mital, Murthy and Chamis 10 used a 
slicing approach to compute the effective elastic constants and microstresses (fiber and matrix stresses) in 
ceramic matrix composites. In this work, a mechanics of materials approach was used to compute the 
effective elastic constants and microstresses in each slice of a unit cell. Laminate theory was then applied 
to obtain the effective elastic constants for the unit cell. Laminate theory was also used to compute the 
effective stresses in each slice, which was used to compute the microstresses. However, in this model 
Poisson effects were neglected in the micromechanics, and the analyses assumed that the material had an 
elastic, strain rate independent deformation. Goldberg extended this slicing approach in order to include 
material nonlinearity and strain rate dependency in the deformation analysis of carbon fiber reinforced 
polymer matrix composites. 1112 This advanced micromechanics model was implemented into classical 
laminated plate theory for the analysis of symmetric thin laminated plates subject to in-plane loading. 13 
In this work, several key issues which are important in the mechanical analysis of composite laminated 
plates under impact loading were considered. However, transverse shear stresses and strains, which are 
also important in the analysis of impact problems, were neglected. 

In the present paper, first the inelastic constitutive model used to model the nonlinear, strain rate 
dependent deformation of the polymer matrix constituent will be briefly described. Next, the previously 
developed micromechanics model 12 will be extended to incorporate the effects of transverse shear stresses 
and strains into the strain rate dependent inelastic composite analyses. Various uniform stress and uniform 
strain assumptions to be applied to the unit cell will be investigated, and the resulting formulation will be 
validated through comparisons to closed form analytical and finite element results. Revisions to the 
micromechanics method to allow for the computation of through the thickness strain and transverse 
Poisson’s ratios will be presented. Finally, the implementation of the revised micromechanics model into 
a refined HOT, in order to improve the ability of the methodology to analyze polymer matrix composites 
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subjected to impact loadings, will be presented and some sample numerical calculations on a 
representative polymer matrix composite will be discussed. 


Constitutive Equations to Analyze Nonlinear Deformation 
of Polymer Matrix Constituent 

To analyze the nonlinear, strain rate dependent deformation of the polymer matrix constituent, the 
Bodner-Partom viscoplastic state variable model, 14 which was originally developed to analyze the 
viscoplastic deformation of metals above one -half of the melting temperature, has been modified. 1 1 In 
state variable models, a single unified strain variable is defined to represent all inelastic strains. 15 
Furthermore, in the state variable approach there is no defined yield stress. Inelastic strains are assumed to 
be present at all values of stress, the inelastic strains are just assumed to be very small in the “elastic” 
range of deformation. State variables, which evolve with stress and inelastic strain, are defined to 
represent the average effects of the deformation mechanisms. 

In the modified Bodner model used in this study, the components of the inelastic strain rate 

tensor, £ - , are defined as a function of the deviatoric stress components S- , the second invariant of the 

deviatoric stress tensor J 2 and an isotropic state variable Z, which represents the resistance to molecular 
flow. The components of the inelastic strain rate are defined as follows 
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where D 0 and n are both material constants, with D 0 representing the maximum inelastic strain rate and n 
controlling the rate dependence of the material. The effective stress, c,„ is defined as 


a e=^2+^ aa kk ( 2 ) 

where a is a state variable controlling the level of the hydrostatic stress effects and c kk is the summation 
of the normal stress components which equals three times the mean stress. The inelastic strains are added 
to the elastic strain tensor to obtain the total strains. 

The rate of evolution of the internal stress state variable Z and the hydrostatic stress effect state 
variable a are defined by the equations 
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where q is a material constant representing the “hardening” rate, and Z x and cq are material constants 
representing the maximum value of Z and a, respectively. The initial values of Z and a are defined by the 

material constants Z 0 and a 0 . The term e\ in equations (3) and (4) represents the effective deviatoric 
inelastic strain rate, which is defined as 
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where ejj are the components of the inelastic strain rate tensor and is the mean inelastic strain rate. In 

many state variable constitutive models developed to analyze the behavior of metals, 15 the total inelastic 
strain and strain rate are used in the evolution laws and are assumed to be equal to their deviatoric values. 
As discussed by Li and Pan, 16 since hydrostatic stresses contribute to the inelastic strains in polymers, 
indicating volumetric effects are present, the mean inelastic strain rate cannot be assumed to be zero, as is 
the case in the inelastic analysis of metals. Further information on the constitutive model, along with the 
procedures required to obtain the material constants, can be found in Goldberg, et al." 


Modification of Micromechanics Model to Include 
Effects of Transverse Shear Deformation 

To compute the effective strain rate dependent, nonlinear deformation response of polymer matrix 
composites based on the response of the individual constituents, a micromechanical model has been 
developed. 12 In this model, the unit cell is defined to consist of a single fiber and its surrounding matrix. 
The composite laminates are assumed to have a periodic, square fiber packing and a perfect interfacial 
bond. Due to symmetry, only one -quarter of the representative unit cell is analyzed. The fiber is assumed 
to be transversely isotropic, linear elastic and rate independent (common assumptions for carbon fibers) 
with a circular cross section. The matrix is assumed to be isotropic with a rate dependent, nonlinear 
deformation response computed using the equations described in the previous section. The unit cell is 
divided into several rectangular, horizontal slices of equal thickness as shown in figure 1 . Due to 
symmetry, only one-quarter of the unit cell needs to be analyzed. All slices are assumed to be in a plane 
stress state, i.e., only in-plane strains and stresses are considered. Each slice then is separated into two 
subslices, one composed of fiber material and the other composed of matrix material. Along the fiber 
direction (direction 1), the strains are assumed to be equal in each subslice of the slice, and the stresses 
are combined using volume averaging. The in-plane transverse normal stresses (direction 2) and in-plane 
shear stresses (direction 12) are assumed to be equal in each subslice of the slice, and the strains are 
combined using volume averaging. Likewise, all the average strains in each slice are assumed to be equal, 
and the stresses are combined using volume averaging. The out-of-plane normal strains (direction 3) are 
assumed to be uniform in each slice, and the volume average of the out-of-plane stresses in each subslice 
is assumed to be zero. The effective properties, stresses and strains for each slice are computed by 
applying these assumptions along with the constitutive equations for the fiber and matrix and solving the 
resulting series of equations. The equivalent stresses, strains and properties for the unit cell can then be 
computed by applying laminate theory to each slice. The advantage of this type of modeling approach 
over other micromechanics methods is that the behavior of each slice is decoupled, and the response of 
each slice can be determined independently. This process results in a series of small matrix equations, 
instead of one large system of equations, which reduces the complexity of the analysis significantly and 
increases the computational efficiency. For elastic materials, the property of transverse isotropy of the 
material could be maintained by using a combination of horizontal and vertical slices and the concept of 
superposition. 9 However, for the nonlinear materials used in this study, superposition cannot be applied, 
and due to the fact that only horizontal slicing is used, the inherent transverse isotropy of the composite 
unit cell is lost. Furthermore, transverse shear effects were neglected in the original analysis. The goal of 
the present work is to include the effects of transverse shear deformation in the micromechanics, taking 
into account the fact that the basic sliced unit cell is not naturally transversely isotropic. 
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Constitutive Models for Fiber and Matrix for Transverse Shear Analysis 

In the development of the revised micromechanics model in which transverse shear effects are 
included, once again the fibers are assumed to be linear elastic and rate independent, and the nonlinear, 
rate dependent behavior of the polymer matrix is computed using the constitutive equations described 
earlier. The constitutive relationships of transverse shear deformation (direction 13 and 23) for the fiber 
and matrix are expressed as 


y F = ( l l G /) aF 6 

T«=(l/Gjo" + y"' 

where "/ and y v/ represent the transverse shear strains of the fiber and matrix, respectively, y Mr represents 
the inelastic transverse shear strain in the matrix, and a F and G W represent the transverse shear stresses of 
the fiber and matrix, respectively. Gy and G m represent the transverse shear modulus of the fiber and 
matrix, respectively. Note that the general form of the equations is the same whether the 13 direction or 
23 direction transverse shear stresses and strains are considered. The only difference would be in which 
shear modulus (13 or 23) for the fiber is used in the equations. 


Different Unit Cell Assumptions for Computing Transverse Shear Deformation 

In this paper, four possible assumptions are chosen for the unit cell in order to compute the local and 
effective transverse shear deformation of the composite ply. To achieve an acceptable level of accuracy in 
the calculations while maintaining a reasonable level of computational efficiency, several relatively 
simple sets of assumptions are proposed and investigated in this work. The various sets of assumptions 
are investigated to determine which result in the most accurate answers. The basic idea of all of these 
assumptions is to assume either a uniform stress distribution or a uniform strain distribution within each 
slice and either the average strain or average stress of each slice is equal in the unit cell. For the slices that 
are composed of both fiber and matrix, two possible assumptions of transverse shear deformation can be 
made in each slice. These are as follows. 

1 . The strains are constant within the slice i and the stresses within the slice are combined using 
volume averaging, that is 


t 

a 1 
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= Y =Y 
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where Vj- is the local fiber volume fraction within slice i. Y and o' represent the equivalent transverse 

strain and stress of slice i, respectively, y 7 and o ' 7 represent the equivalent transverse strain and stress in 
the fiber within slice i, and and o represent the equivalent transverse strain and stress in the matrix 
within slice i. Substituting eq. (6) into eq. (7), the following expression of the equivalent transverse stress 
in slice i in terms of the equivalent transverse strain in slice i and the transverse inelastic strain of matrix 
is obtained 


c'=( v fGf +( l ~ v f )g„ )y' - (i - v} )c m i MI 

where once again the equation is valid for both 13 and 23 direction stresses and strains. 


( 8 ) 
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2. The stresses are constant within the slice i and the strains are combined using volume averaging, 
that is 


a 1 

i 


= a ,F = a lM 

= V' f Y F + (1 - Vj. )y iM 


(9) 


From eqs. (6) and (9), the following expression for the equivalent transverse strain of slice i in terms of 
equivalent transverse stress for slice i and the inelastic strain of the matrix is obtained 

Y = (v‘ /G f +{\ l-V‘ f )/G m )o' +(l-vj \ iMI (10) 

Similarly, two possible assumptions on the transverse shear deformation between the slices of the unit 
cell are investigated. 

1 . The strains are constant between slices and the stresses are combined using volume averaging, 
that is 


y = Y 

o = 2>‘o*. i = l---Nj-+\ 
i 


(ii) 


In eq. (1 1), y and a represent the equivalent transverse shear strain and equivalent transverse shear 
stress of the unit cell, respectively. IV/ represents the number of fiber slices in the quarter of the unit cell 

which is analyzed, and hK represents the ratio of the thickness of slice i to the total thickness of the unit 
cell. 12 

1 . The stresses are constant between slices and the strains are combined using volume averaging, 
that is 


a = a 


Y = X /? fY ! - / = 1- 


•^V / +1 


( 12 ) 


Equations (8), (10), (11) and (12) result in four possible combinations of assumptions to describe the 
transverse shear deformation of the composite unit cell. Table 1 presents the equivalent transverse shear 
modulus G and equivalent inelastic transverse shear strains / of the unit cell resulting by applying the 
various assumptions. These assumptions are denoted by assumption (1) to (4) in this table in order to 
uniquely identify them. 


Numerical Investigation of Effective Transverse Shear Modulus Resulting From Assumptions 

As shown in table 1 , each of the various combinations of assumptions both within the slice and 
between slices can result in significantly different effective transverse shear moduli. To investigate 
qualitatively the differences in effective shear moduli that result by applying the different assumptions, a 
series of numerical computations are performed. In the next section, specific quantitative verification 
studies will be carried out to determine which set of assumptions most accurately predicts the composite 
transverse shear modulus, and by extension the transverse shear response of the composite. 
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The first qualitative analysis is shown in Table 2. In this table, the values of the nondimensionalized 
transverse shear modulus of the unit cell computed using the different assumptions is shown. The 
parameters used in this computation are the nondimensionalized shear modulus of fiber, ( Gf/G m ), which is 
set equal to 15.92 (a typical value for carbon fiber reinforced polymer matrix composites), and the fiber 
volume fraction, Vf, which is equal to 0.6 (again, a typical value). As can be seen in the table, the 
assumptions which are utilized result in significant differences in the effective shear modulus which is 
predicted. Assumption (1) (refer to table 1) results in the highest effective shear modulus by a significant 
amount. The remaining assumptions display some variation in values, with assumption (4) yielding the 
lowest value. The dependence of the transverse shear modulus on founder the various assumptions is 
shown in figure 2. Figure 2 shows the variation of the nondimensionalized transverse shear modulus of 
the unit cell, (G/G m ) with founder the different assumptions. Note that the maximum V f is 0.785 (a 
maximum possible fiber volume fraction for a square pack fiber architecture) and that these figures are 
plotted with a value of ( G f /G m ) = 15.92. The various assumptions used are denoted (1) to (4) as indicated 
in table 1. As seen from figure 2, assumption (1) consistently yields the stiffest modulus, and the variation 
of modulus with fiber volume fraction is linear. The other assumptions yield a nonlinear variation of 
effective shear modulus with fiber volume fraction. Assumption (4) consistently yields the softest 
modulus, and the curves generated by using assumptions (2) and (3) fall in between. The maximum 
deviation between assumption (1) and the remaining assumptions occurs at Vf- 0.6. However, at higher 
volume fractions the differences between the results generated by assumptions (2)-(4) become more 
pronounced, with the maximum difference occurring at Vf = 0.785. Figure 3 presents the variation in 
(G/G,„) with ( Gj!G m ) under the various assumptions, for a constant volume fraction V f of 0.6. From this 
figure, it appears that, for assumptions (2) to (4), the change in effective shear modulus as a function of 
fiber modulus is not very significant, particularly for higher values of the fiber shear modulus. 

Assumption (1) yields a linear increase with increasing fiber modulus. To provide a more general 
representation of the results discussed above, figures 2 and 3 are plotted on the surface of (G/G,„) ~ 
\(Gf/G m ), V/], as shown in figure 4. From these figures, it can be seen that the various assumptions on the 
composite unit cell yield significant differences in the computed effective transverse shear modulus. 


Verification of Assumptions of Transverse Shear Deformation 

As shown in the previous section, by applying different uniform stress and uniform strain 
assumptions to the composite unit cell, both within slices and between slices, significant variations in the 
computed effective transverse shear modulus result. Therefore, it is important to determine which 
assumption yields the most accurate answers in computing the transverse deformations in the 13 and 23 
directions. As discussed earlier, this task is complicated by the fact that the slicing methodology 
employed in this work eliminates the inherent transverse isotropy present in the unit cell. Therefore, the 
appropriate assumptions to use in the transverse shear analysis are determined by comparing the 
equivalent shear modulus computed using the micromechanics method described here (including the 
various assumptions) and modulus values computed using alternative theoretical and numerical methods. 
There are some existing models which can be used to verify these assumptions. A very well known 
model, derived by Hill 17 and Hashin, 18 is used to verify and examine the transverse shear assumptions in 
the 13 direction. To verify and examine the transverse shear assumptions in the 23 direction, the 
commercial finite element software ABAQUS 19 is used to compute the transverse shear modulus. 


Model for Computing Transverse Shear Modulus in 23 Direction in ABAQUS 

Considering the transverse isotopic properties of the unit cell, the model for conducting the transverse 
shear modulus analysis in ABAQUS 19 is reduced to a two-dimensional model. In this procedure, one or 
multiple unit cells are used to calculate the equivalent transverse shear modulus of the unit cell. The 
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elastic constants of the fiber and matrix used in these analyses are shown in table 3, which are 
representative properties for a carbon fiber reinforced polymer matrix composite system and are taken 
from Goldberg, et al." 12 The adequacy of the mesh size is verified via comparison of the results obtained 
by varying the mesh density. When Vf< 0.6, one unit cell is divided into 1600 plane strain elements; and 
when Vf= 0.7, 10000 plane strain elements are used in the calculations. 

When calculating the transverse shear modulus of the unit cell, it is crucial to apply appropriate 
boundary conditions to account for the influence of surrounding unit cells. Determining the boundary 
conditions prior to computation is quite difficult because the edges of unit cell do not remain straight 
during shear deformation. In this paper, an approximate approach is developed by using “stress boundary 
loading” and “strain boundary loading” conditions as shown in figure 5. In stress boundary loading, the 
edges of the unit cell are subjected to tangential displacement loadings (analogous to adding stress 
loading) and the displacements normal to the edges are not constrained. It is important to note that this 
definition does not imply maintaining constant stresses. When strain boundary loading is applied, the 
edges of the unit cell are subjected to combinations of tangential and normal displacement loadings so 
that the deformed shape of the model is a parallelogram. Once again, this does not imply maintaining 
constant strain. In computing the transverse shear modulus, firstly, the stress and strain are calculated at 
the centroid of each element. Then these stresses and strains are combined using volume averaging to 
obtain the average stress and strain, respectively. The effective transverse shear modulus G 23 is calculated 
by using the following expression. 



i / i 


where the summation is over the elements i, V) is the volume of element i, o 23i is the transverse shear 
stress in element i and s 2 3 i is the transverse shear strain in element i. 

Figure 6 illustrates the procedure used to obtain an accurate value of G 22 . At first, by applying stress 

boundary loading on one unit cell, the transverse shear modulus Gi,g is obtained. Then by applying stress 

boundary loading on an assemblage of nine unit cells arranged in a three by three pattern, G 23 calculated 

from the central unit cell is obtained. By applying stress boundary loadings on a larger assemblages of 
unit cells (five by five, seven by seven, etc.), a series of effective transverse shear moduli, 

Goj , G93 , G23 , • • • G23 , are obtained (as shown in fig. 6). As the number of unit cells is increased, the 

computed transverse shear modulus increases. Similarly, by applying strain boundary loadings to a series 
of increasingly larger assemblages of unit cells, a series of effective transverse shear moduli, 

G^3 , G93 , Gjg , • • • G23 , which decrease in magnitude as more unit cells are used, are obtained. When n is 

large enough, while not explicitly shown in Figure 6, G33 ~ G"| ~ G 93 , resulting in a converged 

effective transverse shear modulus of the unit cell. 

Numerical results obtained using this procedure is presented in figure 7. The composite material S2 
Glass/Epoxy, with fiber volume fraction 0.6, is used here. This material was chosen since the glass fiber 
is isotropic, which simplifies the analysis. For this material, both the fiber and matrix materials in this 
composite are isotropic. The shear moduli of Glass and Epoxy are 35.04 GPa and 1.53 GPa, respectively. 
Results in figure 7 show that convergence is obtained using an assemblage of three by three unit cells. It 

can be seen that by applying stress boundary conditions to a unit cell, a value of G^ almost equal to the 

converged multi unit-cell model, G 22 , is obtained. The error is on the order of 0.5 percent, which is very 
small compared to the differences in shear modulus values resulting from the different assumptions. Due 
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to this small difference (between G^ and G 2 3 ), the stress controlled unit cell value of is used to 
conduct the analysis described in the following section. 

Model for Computing Transverse Shear Modulus in 13 Direction 

A solution based on a cylinders model, assuming that the fiber phase is composed of infinitely long 
circular cylinders embedded in a continuous matrix phase, was developed by Hill 17 and Hashin. 18 An 
exact solution of the displacement field was obtained and the equivalent transverse shear modulus in the 
13 direction for the fiber reinforced material was expressed as follows. 

G i3 = G m \G nf (\ + V f ) + G m (1 - V f )]/[G 12 f (1 ■ - V f ) + G m (1 + V f )] (14) 

where G 12 f and G,„ are the shear modulus of the fiber and matrix, respectively, and V f is the fiber volume 
fraction of the unit cell. It must be noted that the unit cell model developed in this work is different from 
that due to Hill 17 and Hashin. 18 However, from a macroscopic point of view, both models result in the 
same equivalent transverse shear modulus for a laminated composite. Thus the results obtained using 
eq. (14) can be used to verify the assumptions discussed earlier. The quantitative analyses and 
determination of the appropriate assumptions to use will be discussed in a later section of the article. 


Study on Through the Thickness Strains and Transverse Poisson’s Ratio 

To accurately compute the through thickness behavior of polymer matrix composites under impact 
loading, modeling through the thickness strains and the transverse Poisson’s ratio become important. In 
the current work, the micromechanics model developed by Goldberg, 12 which did not address these 
issues, is further extended to compute these values. 

Based on the original work reported by Goldberg, 12 each slice in figure 1 is analyzed independently. 
Most of the slices are assumed to have two subslices, with the exception of the top slice. One subslice 
represents fiber material and the other subslice represents matrix material. The top slice is assumed to 
comprise matrix material only. The micromechanics equations presented here are for those slices 
composed of both fiber and matrix material. The stresses in the slices comprising pure matrix can be 
computed using the matrix elastic properties and inelastic constitutive equations (eqs.(l) to (5)) directly. 
The transversely isotropic and isotropic compliance matrices are used to relate the local strains to the 
local stress in the fiber and matrix using the following relations. 
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(15) 


(16) 


where the stresses, Gy, and strains, £y, are defined in a Cartesian frame of reference and Sy represent terms 
in the material compliance matrix. A superscript “7” is used to denote inelastic strains and superscripts 
“F' and “Af ’ are used to indicate fiber subslice and matrix subslice, respectively. Subscripts “f' and “m” 
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refer to fiber and matrix related properties, respectively. Out-of-plane normal stresses are included in the 
analysis since although each slice is assumed to be in a state of global plane stress, the individual 
subslices are assumed to be in a full three-dimensional state of stress. The relationships between shear 
stresses and shear strains are not presented here because the through the thickness strain and transverse 
Poisson’s ratio are independent of these quantities. 

Based on the slice geometry (fig. 1), slice fiber volume fraction and the basic definitions governing 
displacement and force equilibrium, the uniform stress and uniform strain assumptions used for the 
micromechanics model can be stated as 


£ 22 = V f e 22f +( 1 -V})e , 22 m 
e 33/ = £ 33 m =£ 33 

... . . (17) 

a ll = V f c llf +(!-! / /)°h m 
C 22 / = a 22 m = a 22 

a' 3 =0 = V i f o i 33f +(\-V i f )o i 33m 

where the superscript i represents variables related to slice i. 

Combining the uniform stress and uniform strain assumptions (eq. (17)) with the constituent stress- 

strain relations (eqs. (15) and (16)), the stress variables, c'\ \ , , o' 33 , and can be derived in 

terms of the total and inelastic strains in the 1 and 2 directions as follows 


,-iF _ , > iA „ i i / ) iA „ i , c iAI 
U 11 _ ^ ^12 fc 22 

r-iM _ (\iKpi i /I'/tp' j.p® 

U 11 _ ^ll fc ll _r ^12 fc 22 " r fc l 

C 22 = ^2i e i 1 + ^22 £ 22 +£ 2 

_ iF _ , > iA „ i , / > iA „ / , _ iAI 
u 33 _ ^31 fc ll ^32 fc 22 _r fc 3 



= QmZ 1 U +Q\ , 7 e ?9 +£ 


iB~i 


JBI 


'3111 


32 22 


(18) 


where Q' are the effective stiffness matrix terms and the various e ,/ are collections of the inelastic strain 
terms. 12 Combining eqs. (17) and (18), the following expression for the out-of-plane normal strain in slice 
i can be expressed (in terms of the computed coefficients H') as follows. 


'33 


= H[ l e u +e 


12 c 22 


il 


(19) 


where e 3 represents collected inelastic strain terms of slice i. Similar expressions can be derived for the 

slices containing matrix only. Since the normal stresses in the 3 direction are assumed to be constant 
(zero) through the thickness direction, the through the thickness normal strain, e 33 , is obtained by 
combining the effective strain values at each slice (eq. (19)) using the following expression. 

N f + 1 

£ 33 - Z ^} £ 33 (20) 

i=l 


NAS A/TM— 2004-2 1 3420 


10 



In eq. (20), h K is the thickness ratio for slice i, which is defined as h 1 / h T , where h' is the thickness 

of slice i, h T is the total thickness of the quarter of the unit cell which is analyzed and N f is related to the 
number of fiber slices as mentioned before. Substitution of eq. (19) into eq. (20) yields the following. 


-33 


~ H n £ n +H n £ 22 +e 3 


( 21 ) 


where 


N f+ l 

H U = Z h' f H It 
1=1 
N f + 1 

H 12 - X h l f H [ 2 (22) 

i = 1 
A (f+l 

£ 3 = E /7 / £ 3 
j=l 

To determine the effective inelastic strain in the through thickness direction ( e 33 )’ e T ( 21 ) is 
rearranged as 

£ 33 = ^ll( £ ll _£ 11 ) + ^12 ( £ 22 — £ 22 ) + £ 33 ( 22 ) 

Using eqs. (20) and (21), the inelastic normal strain can be expressed as follows. 

£ 33 = ^11 £ 11 + ^12 £ 22 +e 3 ( 24 ) 

To determine the transverse Poisson’s ratios, the procedure proposed by Whitney 9 is used. Equation 
(21) is simplified for elastic loading along two directions. First, for loading along the 1 direction, using 
the usual definition for Poisson’s ratio, v , the following equation is obtained. 

~ V 13 £ 11 = ^11 £ 11 -V 12^12 e ll ( 2 ^) 

which can be simplified to yield the following expression 

v 13 = -//u + v l 2 //t2 (26) 

Likewise, for loading along the 2 direction, the following result is obtained. 

- v 23 £ 22 = H\2 £ 22 _ v 21 ^ 11 £ 22 ( 27 ) 

which can be simplified to yield 

v 23 = ~//i2 + v 21-^ll (28) 
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Implementation of Higher Order Laminated Plate Theory 
Into Strain Rate Dependent Micromechanics Model 

A micro-macro analysis procedure is developed by coupling the strain rate dependent 
micromechanical model with a refined HOT. In this procedure, the verified assumptions which appear to 
provide the best descriptions of the transverse shear deformation, which are discussed in a later section of 
this article, are used. On the laminate level, the traditional HOT 3 4 is refined to capture the inelastic strain 
effects. In this work, hydrothermal effects are neglected. 


Modified Higher Order Laminated Plate Theory 

In the higher order theory, cubic through-the -thickness variations are assumed to describe the inplane 
deformations and the out of plane deformation is assumed to be independent of laminate thickness. This 
allows an accurate and efficient description of transverse shear stresses, which are important in 
anisotropic laminates. The displacement field is expressed as follows. 

u(x, y,z ) = Uq(x, y) + z$ x (x, y) + z 2 <p x {x, y)+z 3 Q x (x, y) 

t>U, y , Z ) = t>0 U> y ) + zty y (x, y ) + Z 2 9 v (x, y ) + r 3 6 v (x, y ) (29) 

w{x,y,z)=w 0 {x,y) 


where u(x,y,z), v(x,y,z), w(x,y,z) are the displacements of a point (x,y,z) and u 0 (x,y), o 0 (x,y), w 0 (x,y) are the 
corresponding values in the mid-plane and z is the distance from the mid-plane along the thickness. 
Application of the stress-free boundary conditions, a 13 | 7= +/,/? = 093 ] , =+/l/9 = 0 , at top and bottom 

surfaces, results in the following expressions for strains associated with the displacements, described in 
eq. (29) . 
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where 
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( 31 ) 


In eqs. (30) and (31), the strains are expressed using conventional engineering notations (that is, 

1 = 11,2 = 22, 6 = 12, 4 = 23 and 5 = 13) and h represents the plate thickness. 

Combining eqs. (30) and (31) and the equivalent constitutive relationship of each ply, the resultant 
forces are expressed in terms of the mid-plane strains and inelastic resultant forces. For a symmetric 
laminate, these expressions can be written as follows 


where 
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(32) 


(33) 


In eqs. (32) and (33), all variables are defined in the structural axis system and conventional 
engineering notations are used to express normal and shear stresses and strains. The quantities IV) and Q, 

are the in-plane force resultants and transverse shear force resultants, respectively and N- and Q- are 
the inelastic in-plane force resultants and inelastic transverse shear force resultants, respectively. Ay and 
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Djj are laminate stiffness matrix components and Q- , a, and e y are the stiffness matrix components, 

stresses and the inelastic strains of each ply in the structural axis system, respectively. In this refined 
HOT, both transverse shear effects and inelastic effects which are important in the impact problem are 
considered. The in-plane force resultants and transverse shear force resultants are decoupled in a 
symmetric laminate. It must be noted that the Q used here have a different meaning from those used in 
eq. (18) which represent the calculated effective stiffness matrix terms of slice i in a unit cell. 


Computational Procedure 

The details of the computational procedure are illustrated in figure 8, where the left half represents the 
macro-component of the analysis procedure using the higher order theory (HOT), and the right half 
describes the micro-component of the analysis based on the extended micromechanics model. First, the 
mid-plane strains of the laminated plate are obtained from direct input or a finite element program. Then 
the strains of each layer are calculated using eq. (30). Because each layer is represented by a unit cell, the 
strains of each layer are transferred into a unit cell in the micro-part of this procedure. In the unit cell, the 
input strains are used to calculate the stresses of each subslice. Then the inelastic strains of the matrix 
subslice are calculated using eqs. (1) to (5). This is followed by calculation of the equivalent inelastic 
strains of the unit cell, which are then transferred back to the macro-part of this procedure. Finally the 
resultant forces of the laminate are calculated using mid-plane strains of the laminate and the inelastic 
strains of each layer. 


Results 

First, numerical analysis is conducted to investigate the validity and accuracy of the different 
assumptions discussed in the Modification of Micromechanics Model to Include Effects of Transverse 
Shear Deformation section (table 1). A representative polymer matrix composite material, IM7/977-2, 
with a fiber volume fraction of 0.6, is used in this study. Detailed convergence studies on Nf was 
conducted in earlier work. 12 Therefore, in present work, 5 fiber slices (Nf- 5) in the quarter of the unit 
cell are used. The uniform stress and uniform strain assumptions that yield the best values to describe the 
deformations in the 13 and 23 directions are determined and these are subsequently used in the structural 
analysis of composite laminated plates. The mechanical behavior of composite plates with various 
stacking sequence are investigated under various strain rate loadings. 


Verification of the Transverse Shear Modulus G23 and G13 

For the verification studies, properties for the carbon fiber reinforced polymer matrix composite 
system, IM7/977-2, presented by Goldberg, et al. 11 ' 12 were used. The material properties are listed in 
table 3. Only the elastic transverse shear properties of the fiber and matrix are used here. 

Comparison of the effective transverse shear modulus values, G23, obtained using the four 
assumptions on the unit cell discussed earlier (and shown in table 1) with those obtained using ABAQUS 
are presented in figure 9. Figure 9 represents the variation in transverse shear modulus G23 with V f 
(overall composite fiber volume fraction). From this figure, it can be concluded that the results obtained 
using assumption 3 have a good correlation with those obtained using ABAQUS when the fiber volume 
fraction is less than 0.3. When V f is equal to or greater than 0.6, the shear modulus results obtained using 
assumption 4 show better correlation with ABAQUS. Since the normal range of V/for composite 
materials used in impact applications is typically between 0.5 and 0.7, assumption 4 is chosen in this 
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work. The transverse shear deformation in the 23 direction is obtained by applying this assumption to a 
unit cell and these values are used in subsequent computations. 

The effective transverse shear modulus, Gn, obtained using the four assumptions and the values 
obtained using the formulation from Hill 17 and Hashin 18 (eq.14) are compared in figure 10. As can be seen 
in the figure, the results obtained using assumption 3 show a good correlation with the results obtained 
using the closed form formulation. At high volume fractions, due to differing assumptions as to how the 
unit cell is constructed, there is some deviation between the analytical result and the values computed 
using the micromechanics model. Since the fiber volume fractions of greatest interest for impact 
applications range between 0.5 and 0.7, assumption 3 is selected to describe the transverse shear 
deformation in the 13 direction. 


Mechanical Response of Various Laminate Plates Under Different Strain Rate Loadings 

In this work, IM7/977-2 composite laminated plates with various stacking sequences are used to 
investigate the material response under various loading directions and various strain rate loadings. The 
thicknesses of all the laminates considered are 0.01 m. Material properties of the fiber and matrix are 
shown in table 3. The fiber is assumed to be linear elastic and the matrix is modeled using the rate 
dependent, nonlinear model and material properties summarized earlier in this article and described in 
detail in Goldberg, et al. 11 The composite material has a fiber volume fraction of 0.6. Five fiber slices are 
used (Nf= 5) in the quarter of the unit cell which is analyzed. Assumptions (3) and (4), which appear to 
provide reasonable descriptions of the transverse shear deformation in the 13 and 23 directions, 
respectively, are used. Here, directions 1 , 2 and 3 represent the material coordinate system and directions 

x, y and z represent the laminate coordinate system. Subscripts x and y represent in-plane directions and z 
represents the out-of-plane direction. Thus, for the unidirectional composite plate (stacking sequence [0]), 
the material axis system coincides with the structural axis system. 

The mechanical response of a composite plate with [0]4 stacking sequence is shown in figure 11. This 
figure represents different stress-strain curves under different strain rate uniaxial loadings in the x, y, xz 
and yz directions. The experimental values (scatter points) in figure 1 lb and 12a were taken from 
Goldberg, et al. 1 1 The high strain rate testing was conducted using a Hopkinson bar and the low strain rate 
testing was conducted using an Instron machine. In figure 1 1(b), the apparent nonlinearity in the 
experimental high strain rate stress-strain curve along the y direction, and the difference between the 
experimental and computed results at the high strain rate, is most likely due to the lack of stress 
equilibrium in the experiments at low strains and the specimen geometry used for the experimental tests. 
More details regarding the experimental procedures and results are provided in Goldberg et al. 1 1 
Unfortunately, the experimental values were limited to in-plane directions only. The significance of strain 
rate effects is investigated by studying the effects on the transverse shear stress-strain curves. More 
experiments along the transverse directions will be conducted in the future. Figure 1 1(b) and figure 12(a) 
(discussed later) show good correlations between the experimental values and the computational results 
under various stacking sequences and various strain rate loadings. Although, the comparisons are based 
on limited available experimental values, some meaningful observations can still be made regarding the 
accuracy of the modeling technique. Based on these comparisons, it appears that the developed micro- 
macro mechanical framework can accurately describe the structural behavior of rate dependent composite 
laminated plates under various strain rate loadings. From figure 1 1, it can be seen that almost no strain 
rate effect is exhibited along the fiber direction, as expected. However significant rate dependencies in the 

y, xz and yz directions are observed. The stresses increase with an increase in strain rate. Physically, this 
is due to the fact that in these directions the response of the matrix plays a more dominant role in the 
stress analysis. In case of the unidirectional laminate, since the material axis and the structural axis 
coincide, this material rate dependent behavior is reflected in the laminate stresses (in y, xz and yz 
directions). An important observation is the fact that the magnitudes of the transverse shear stresses are 
nontrivial. The magnitude of the transverse stresses will also increase with an increase in laminate 
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thickness. Thus, nontrivial, nonlinear and rate dependent transverse shear stresses, which are important 
for an impact problem, are observed by using the new modeling techniques. Therefore, it is expected that 
the developed refined micro-macro mechanical model, which includes transverse shear deformation, will 
significantly improve the level of accuracy in the impact analysis of composite laminates. 

The mechanical response of a composite plate with a [45] 4 stacking sequence under uniaxial loadings 
along x and xz directions is shown in figure 12. Significant rate dependency can be observed both in the 
normal and in the transverse shear directions. This is because, in case of an orthotropic laminate, the 
applied loading can be decomposed into loading along the fiber direction and loading along the in-plane 
transverse normal and shear directions. The loading along fiber direction will not exhibit strain-rate 
dependency. However, the stresses in the in-plane transverse normal and shear directions will exhibit 
significant differences under different strain rate loadings. This combined effect results in the rate 
dependency shown in figure 12 for the [45] 4 laminate. These results indicate that along with the material 
properties of the laminate being considered, plate geometry, stacking sequence and loading directions will 
also influence the rate dependency in the structural response of composite laminated plates. 

Figure 13 shows the mechanical response along xz and yz directions of a composite plate with a 
[90/0] s stacking sequence. It can be seen that the transverse shear responses are different along the two 
directions. This can be explained as follows. In HOT, the assumed inplane displacements (u, v) are cubic 
functions of the thickness coordinate, which results in parabolic distributions of transverse shear strain 
along the thickness direction. Therefore, in transverse shear stress analysis, the layers closer to mid-plane 
play a more important role than those near the top and bottom surfaces. This result in different values of 
‘average’ transverse shear moduli of the laminate along yz and xz directions. As a result, the stress-strain 
curves are different along the two directions. If FSDT is used, due to the linear distribution of 
displacements (u, v), the transverse shear strains remain constant through the thickness. Therefore, the 
locations of the layers do not affect the transverse shear stress analysis. As a result, in the case of a cross- 
ply laminate, the effects shown in figure 13 cannot be observed using FSDT. 


Conclusion 

An investigation on the nonlinear responses of polymer composite laminate plates under different 
strain rate loadings is conducted. A strain rate dependent micromechanics model is further extended to 
account for the transverse shear effects that are important in the impact problem. Four different 
assumptions of transverse shear deformation are investigated to improve the developed strain rate 
dependent micromechanics model. Theoretical and finite element investigations are conducted to validate 
these assumptions. Two of these assumptions are found to be adequate in describing the transverse shear 
deformation along the transverse directions, over the normal range of fiber volume fraction. Based on 
these assumptions, a method to determine through the thickness strain and transverse Poisson's ratio is 
developed. A higher order laminated plate theory, which can capture the through-the-thickness transverse 
shear stresses, is extended to include inelastic effects and is used in the laminate analysis. An integrated 
micro-macro numerical procedure is developed by implementing the revised micromechanics model into 
the refined higher order laminated plate theory, which now contains the inelastic strain effects. The 
mechanical response of IM7/977-2 composite laminate plates with different stacking sequences are 
investigated under various strain rate loading conditions. The results illustrate the significant strain rate 
dependency and nonlinear effects even for moderate strain rate loadings. The strain rate dependency is 
also influenced by laminate geometry and ply stacking sequence. The results show that transverse shear 
stresses play an important role at both macro and micro level behavior of composites. In particular, these 
effects cannot be neglected in the impact analysis of composites laminates. 
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Table 1. — Formulations under different transverse shear assumptions 
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Table 2. — Nondim ensionalized shear modulus of unit cell under diff erent assumptions 


Assumption 

G/G m 

1 

9.95 

2 

3.11 

3 

4.78 

4 

2.31 


Table 3. — Material properties of the fiber and matrix 


IM7 fibers 

Ell(Gpa) 

E22(Gpa) 

G12(Gpa) 

G23(Gpa) 

vl2 

v23 




276 

13.8 

20 

12.4 

0.25 

0.25 




977-2 polymer matrix 

e (/sec) 

E(Gpa) 

V 

Do(/sec) 

n 

Zo(Mpa) 

Zj/Mpa) 

q 
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ai 

9E-5 

3.52 

0.4 

1.00E+06 

0.8515 

259.496 

1131.371 
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0.1289 
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1.9 

3.52 
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6.33 
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Figure 1 . Schematic showing relationship 
between unit cell and slices. 



Figure 2. The variation of G/Gm with Vf. 



Figure 3. — The variation of G/Gm with Gf/GM 
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Figure 4. — The variation of G/Gm with 
Gf/GM and Vf. 
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Figure 5. — Stress and strain boundary condition. 


Figure 6. — Calculation of transverse 
shear modulus, G23. 
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Figure 7. — Shear moduli under different 
boundary loadings. 


Figure 8. — Illustration of analysis procedure. 
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Figure 9. — Variation of G23 with Vf. Figure 10. — Variation of G13 with Vf. 
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Figure 11. — Comparison of stress-strain curves of composite with stacking sequence [0°]2s at different strain 
rate loadings, (a) In x direction; (b) In y direction; (c) In xz direction; (d) In yz direction. 
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Figure 12. — Comparison of stress-strain curves of composite with stacking sequence [45°]2s at different strain 
rate loadings, (a) In x direction; (b) In xz direction. 
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Figure 13. — Comparison of stress-strain curves of composite with stacking sequence [0°/90°]s at different 
strain rate loadings, (a) hr xz direction; (b) hi yz direction. 
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